﻿#include <algorithm>
#include <gtest/gtest.h>
#include <iostream>
#include <vector>

#include <xy/alg/interpolator.h>

using namespace xy;
using namespace xy::geom;
using namespace xy::alg;
using namespace xy::math;

TEST(TestInterpolator, Newton)
{
    // 插值算法
    {
        std::mt19937 generator((unsigned int)time(nullptr));      // 使用 time(nullptr) 作为随机数生成器的种子
        std::uniform_real_distribution<double> uniform(0, 100.0); // 在 [0,1] 区间内的均匀分布

        std::vector<double> x, y;

        int N = 10;
        for (int i = 0; i < N; i++)
        {
            x.push_back(uniform(generator));
            y.push_back(uniform(generator));
        }

        // 测试插值函数和 NA 算法
        interpolator_newton<double> interp(x, y);
        auto f = interp.interpolate();
        for (int i = 0; i < N; i++)
        {
            // 计算相对误差，拉格朗日插值不稳定
            EXPECT_LE(std::abs(f.value(x[i]) - y[i]) / std::abs(y[i]), 1e-3);
            EXPECT_LE(std::abs(interp.value(x[i]) - y[i]) / std::abs(y[i]), 1e-3);
            EXPECT_FLOAT_EQ(interp.neville_aitken(x, y, x[i]), y[i]);
        }
    }
}

TEST(TestInterpolator, CubicBSpline)
{
    auto f = [](double x) -> vec3 { return vec3{8.0, x, std::abs(x) + x / 2 - x * x}; };

    std::vector<double> x;
    std::vector<vec3> y;

    int N = 10;
    for (int i = 0; i < N; i++)
    {
        x.push_back(-1 + 2 * i / (N - 1.0));
        y.push_back(f(x[i]));
    }

    interpolator_cubic_bspline<3> bspline;

    {
        auto complete =
            bspline.interpolate(y, bspline.Complete, bspline.Complete, vec3{1.0, 0.0, 0.0}, vec3{1.0, 1.0, 1.0});
        auto [first, last] = complete->bound();
        auto [p1, d1] = complete->d1(first);
        auto [p2, d2] = complete->d1(last);
        EXPECT_LE((d1 - vec3{1.0, 0.0, 0.0}).modulus(), numeric::distance);
        EXPECT_LE((d2 - vec3{1.0, 1.0, 1.0}).modulus(), numeric::distance);
    }

    {
        auto natural =
            bspline.interpolate(y, bspline.Natural, bspline.Natural, vec3{1.0, 0.0, 0.0}, vec3{1.0, 1.0, 1.0});
        auto [first, last] = natural->bound();
        auto [p1, d1, dd1] = natural->d2(first);
        auto [p2, d2, dd2] = natural->d2(last);
        EXPECT_LE((dd1 - vec3{1.0, 0.0, 0.0}).modulus(), numeric::distance);
        EXPECT_LE((dd2 - vec3{1.0, 1.0, 1.0}).modulus(), numeric::distance);
    }

    // 曲面插值
    {
        auto f = [](double u, double v) -> vec3 { return vec3{u, v, sin(u + v) * cos(u - v)}; };

        int N = 10;
        int M = 15;
        std::vector<std::vector<vec3>> points(N, std::vector<vec3>(M));
        for (int i = 0; i < N; i++)
            for (int j = 0; j < M; j++)
                points[i][j] = f(5.0 / (N - 1) * i, 5.0 / (M - 1) * j);
        auto surface = interpolator_cubic_bspline<3>::interpolate(points);
    }
}